function coe = State2Coe(S,u)

s = size(S);
%[]Returns the dimensions of the state matrix.

coe = zeros(s);
%[]Allocates memory for the classical orbital elements matrix.

I = [1; 0; 0];
%[]Inertial X-axis vector.

K = [0; 0; 1];
%[]Inertial Z-axis vector.

for k = 1:s(2)
    
    R = S(1:3,k);
    %[km]Current position vector.
    
    V = S(4:6,k);
    %[km/s]Current velocity vector.
    
    H = cross(R,V);
    %[km^2/s]Current specific angular momentum vector.
    
    N = cross(K,H);
    %[km^2/s]Current ascending node vector.
    
    r = norm(R);
    %[km]Magnitude of the current position vector.
    
    v = norm(V);
    %[km/s]Magnitude of the current velocity vector.
    
    Dot = dot(R,V);
    %[km^2/s]Dot product of the current position and velocity vectors.
    
    E = 1 / u * ((v^2 - u / r) * R - Dot * V);
    %[]Current eccentricity vector.
    
    h = norm(H);
    %[km^2/s]Magnitude of the current specific angular momentum.
    
    n = norm(N);
    %[m^2/s]Current ascending node vector magnitude.
    
    epsilon = v^2 / 2 - u / r;
    %[km^2/s^2]Current specific mechanical energy.
    
    coe(1,k) = -u / (2 * epsilon);
    %[km]Current semimajor axis.
    
    coe(2,k) = norm(E);
    %[]Current eccentricity magnitude.
    
    coe(3,k) = acosd(dot(K,H) / h);
    %[deg]Current inclination.
    
    coe(4,k) = acosd(dot(I,N) / n);
    %[deg]Current right ascension of the ascending node.
    
    if N(2) < 0
        
        coe(4,k) = 360 - coe(4,k);
        %[deg]Right ascension of the ascending node quadrant check.
        
    end
    
    coe(5,k) = acosd(dot(N,E) / (n * coe(2,k)));
    %[deg]Current argument of periapsis.
    
    if E(3) < 0
        
        coe(5,k) = 360 - coe(5,k);
        %[deg]Argument of periapsis quadrant check.
        
    end
    
    coe(6,k) = acosd(dot(E,R) / (coe(2,k) * r));
    %[deg]Current true anomaly.
    
    
    if Dot < 0
        
        coe(6,k) = 360 - coe(6,k);
        %[deg]True anomaly angle quadrant check.
        
    end
    
end


coe = real(coe);

end
%===================================================================================================